Velocity correlations in the dense granular shear flows: 
Eff^ects on energy dissipation and normal stress 
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We study the effect of pre-coUisional velocity correlations on granular shear flow by molecular 
dynamics simulations of the inelastic hard sphere system. Comparison of the simulations with the 
kinetic theory reveals that the theory overestimates both the energy dissipation rate and the normal 
stress in the dense flow region. We find that the relative normal velocity of colliding particles is 
smaller than that expected from random collisions, and the discrepancies in the dissipation and 
the normal stress can be adjusted by introducing the idea of the coUisional temperature, from 
which we conclude that the velocity correlation neglected in the kinetic theory is responsible for 
the discrepancies. Our analysis of the distributions of the pre-coUisional velocity suggests that the 
correlation grows through multiple inelastic collisions during the time scale of the inverse of the 
shear rate. As for the shear stress, the discrepancy is also found in the dense region, but it depends 
strongly on the particle inelasticity. 

PACS numbers: 47.57.Gc,45.70.Mg,47.45.Ab,83.10.Rs 
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I. INTRODUCTION 



Granular media can flow like a fluid under a certain 
situation. In the case of the rapid granular flow, where 
the density is relatively low and interactions are domi- 
nated by the instantaneous collisions, the kinetic theory 
of dense gases [l| is extended to the inelastic hard spheres 
to derive the constitutive relations . In the theory, the 
density correlations is taken into account to some extent 
but not the velocity correlations in most of the cases. As 
the flow gets denser, however, the molecular chaos as- 
sumption becomes questionable. In addition, the inter- 
actions may no longer be approximated by the instanta- 
neous collisions but enduring contacts take place around 
the random closed packing fraction. The comprehensive 
granular rheology including the rather complicated dense 
regime has not been established yet. 

During the last several years, careful experiments and 
large-scale molecular dynamics simulations have been 
done on the dense granular flows [1, |3j S @1- One of 
the important model systems that has been intensively 
studied is the steady flow down a slope under the gravity, 
where we can control the ratio of the shear stress S to the 
normal stress N by changing the inclination angle 9. In 
this system, it has been found that the packing fraction 
v in the bulk of the flow is constant and is determined 
solely by the inclination angle 9; in other words, v is in- 
dependent of the total flow hight H and/or the roughness 
of the slope HQ. 

This interesting feature has been qualitatively under- 
stood by using the Bagnold Scaling ^Tj], which states the 
shear stress S is proportional to the square of the shear 



rate 7: 



(1) 



Here, m is the particle mass, and a is the particle di- 
ameter. This scaling can be understood by dimensional 
analysis of the rigid granular flow, where the inverse of 
the shear rate 7"^ is the only time scale in the system. 
This scaling applies to the normal stress N also, which 
gives 



N = 



(2) 



In the slope flow under gravity, the force balance gives 
S/N = tan 6*. Thus we finally have 



N 



B{u) 



~ tan( 



(3) 
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i.e., the packing fraction v is determined by the inclina- 
tion angle 9. 

This dimensional analysis does not hold when the time 
scales other than 7"^ come into the problem, e.g. the 
time scales of the particle deformation ^8], but not only 
the constant density profile but also the Bagnold scal- 
ing itself has been found in the numerical simulations of 
dense steady fiow down a slope for hard enough parti- 
cles |5j]. 

In the slope flow simulations, the value of the packing 
fraction v has been shown to increase upon decreasing 
the inclination angle 9, and eventually the flow stops at 
a finite angle 0stop; namely, A{y'^B{i') is an decreasing 
function of v in the dense region 5,'6'|. One can interpret 
the transition at ^stop as the jamming transition [3, 9]. 

The theoretical analysis of the functional form of 
A{v) / B{v) has been done by Louge [lO] using the ki- 
netic theory, but he found the opposite dependence in the 
dense region, namely, the theory gives increasing packing 
fraction v upon increasing inclination angle 9 as shown in 
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briefly summarize the inelastic hard sphere model and 
the constitutive relations based on the kinetic theory. We 
summarize our simulation method and present the results 
in section Hm The discussion and the summary are given 
in section HVl 



II. 



INELASTIC HARD SPHERE MODEL AND 
THE KINETIC THEORY 



FIG. 1: (color online) The ratio of the shear stress to the 
normal stress S/N versus the packing fraction i/ from the 
simulation data (ep = 0.70 (■), 0.92 (•), and 0.98 (o)) and 
the plot of Eq. (|20p from the kinetic theory (ep = 0.70 (* 
connected by dashed line), 0.92 (+ connected by dashed line), 
and 0.98 (x connected by dashed line). For the simulation 
data, the average normal stress A'' = ^{N^ + Ny + N^) is 
used. 



Fig.[l] where curves from a kinetic theory |ll| is shown by 
symbols connected by dashed lines for the various resti- 
tution coefficients Cp. 

Several explanations for this discrepancy have been 
proposed, such as the enduring contact [13, , the Bur- 
nett order (the second order of the spatial gradients) ef- 
fect and the particle roughness etc., but the 
subject is still under debate. 

Recently, the present authors Q have made a detailed 
comparison between the simulation results of the dense 
slope flow and the kinetic theory by Jenkins and Rich- 
man [l3|. In contrast with the rather good agreement 
for the stresses, it has been found that the kinetic the- 
ory overestimate the energy dissipation rate F, and this 
discrepancy is responsible for the contradicting behavior 
in the kinetic theory, i.e. A{v)/B{v) increases with the 
packing fraction u. 

The authors conjectured that the discrepancy in the 
energy dissipation rate T should be caused by the veloc- 
ity correlations enhanced by the inelastic collisions; the 
decrease of the relative normal velocity through the in- 
elastic collisions results in reduction of the energy loss 
per collision. Such an effect has been noticed in granular 
gas simulations without shear [l^, , and the velocity 
correlations has been investigated analytically [H, [l^ . 

However, the situation is rather complicated under 
shear, because the shear tends to break the correlations. 
The spatial velocity correlations in granular flow under 
shear has not been carefully studied so far [20j . 

In this paper, we study the velocity correlation in the 
sheared granular flow, focusing its effects on the energy 
dissipation rate and the stress. We adopt the simple 
shear flow of the inelastic hard spheres as a model system, 
in accordance with most of the kinetic theory analysis. 
Note that the enduring contacts is not allowed in the hard 
sphere model, whose effects are often under debate in the 
soft-sphere model simulations of the slope flow 0] . 

This paper is organized as follows. In section [TTl we 



The inelastic hard sphere model is one of the simplest 
and widely-used models of granular materials [2, 22] . The 
particles are infinitely rigid, and they interact through 
instantaneous two-body collisions. We adopt the sim- 
plest collision rule for the monodisperse smooth hard 
spheres with diameter cr, mass m, and a constant nor- 
mal restitution coefficient ep in three dimensions as fol- 
lows: The particle i at the position with the veloc- 
ity Ci collides with the particle j \i \ri — rj\ = a and 
(ri — Tj) ■ {ci — Cj) < 0, and their post-coUisional veloci- 
ties c* and c* are given by 



Ci - • {a -Cj)] n,j, (4) 



2 

1 + e 



- ["ij ■ (Ci - Cj)]nij, (5) 



respectively. Here, riij is a unit vector defined as riy = 
— rj)/\ri ~ rj\. The collision is elastic when Cp = 1, 
and inelastic when < Cp < 1. In the inelastic case, the 
particles lose the kinetic energy every time they collide, 
thus external drive is necessary to keep particles flowing. 

We compare the simulation results of the inelastic hard 
spheres with the constitutive relations obtained from the 
Chapman-Enskog method [l| , which has been developed 
in the kinetic theory of gases. In this paper, we employ 
those by Garzo and Dufty pd| . who have improved the 
previous studies [1, [13, ; that is limited to the weakly 
inelastic case ((1 — Cp) <C 1), to include the case with any 
value of the restitution constant Cp under the assump- 
tion that the state is near the local homogeneous cooling 
state 



24\. 



In the following, we briefly summarize the kinetic the- 
ory to derive the constitutive relations. The hydrody- 
namic variables are the number density field n{r,t), the 
velocity field u{r,t), and the granular temperature field 
T{r,t), defined in terms of the single-particle distribution 
function f{r,c,t) as 



(6) 
(7) 



n{r,t) = J f{r,c,t)dc, 
u{r,t) — — J cf{r,c,t)dc, 
T{r,t) = £ J{c-urf{r,c,t)dc. (8) 
The hydrodynamic equations for these variables are given 
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TABLE I: The dimensionless functions in the constitutive re- 
lations from Ref. [Till. 



/i(^) 

C 



Xl + 2(l + ep)^(7oH) 
[^'=*(l + |,.ffoM(l + e,)) + |7*] 



72(1- 



(0)* 



^3/2^ 

ffoM(l-i(l-ep)^)(l-^c*(ep)) 
30^^(1 -e^)(H-Ac*(ep)) 
^!.2goM(l + ep)(l-^c*(ep)) 



c*(ep) 32(1 - ep)(l - 2e^)[81 - 17ep -^ 30e^(l e^)]- 



by 



dn ^ , . 

mn— — h mnu ■ vu 
at 



dT 



VT j 



0, (9) 
-VS, (10) 
-V ■q-T.-.E-r, (11) 



where S is the stress tensor, q is the heat flux, and 
E is the symmetrized velocity gradient tensor: Eij = 
^{■£rUj + -£rUi)- Note that the energy dissipation rate 
r in Eq. ([TT|) appears due to the energy loss through the 
inelastic collisions, which gives peculiar features to the 
granular hydrodynamics. 

The constitutive relations for S, q, and F are deter- 
mined by the single-particle distribution f{r,c,t). Its 
time evolution depends on the two-particle distribution 
function /^^^(ri, r2, Ci, C2, i) through the two-particle 
collision; the n-particle distribution function depends on 
the {n + l)-particle distribution function. This is known 
as the BBGKY hierarchy [H. 

In the Enskog approximation, the two-particle distri- 
bution at collision is approximated as 

/^^^(ri,ri + an2i,Ci,C2,t) 
= 9oMf{ri,Ci,t)f{ri+<Tn2i,C2,t) (12) 

to close the BBGKY hierarchy at the single-particle dis- 
tribution [T], [ll[. Here, 50 (i^) is the radial distribution 
function at distance a, and depends on the packing frac- 
tion v — ^ncr^n. The term 170(1^) represents the positional 
correlations, and the actual procedure to determine the 
functional form of go{i') is presented in subsection llll B l] 
The correlations in the particle velocities are neglected 
under the molecular chaos assumption. 

The constitutive relations for the hydrodynamic equa- 
tions have been obtained in ref. [IH by the Chapman- 
Esnkog method with the approximation (|12p up to the 
Navier-Stokes order (i.e. the first order of the spatial gra- 
dients). In the simple steady shear flow with constant n, 
constant T, and u{r) = (72,0,0), the nonzero terms are 
the pressure, or the normal stress 



the shear stress 

S = E,,, = S,,, = Siiy,T) = m^/^a-^f2iiy)VTj, (14) 

and the energy dissipation rate 

r - T{iy, T) = m-^'^a-^h{v)T^'^. (15) 

The dimensionless functions fi{v) are listed in Table U 

In the simple shear flow, Eqs. ^ and (jTU]) are auto- 
matically satisfied with the constant normal stress N and 
the constant shear stress S. The energy balance equation 
PT|) gives 



S'7 - r = 0, 



(16) 



because there is no heat flux q. Equation (fT6|l means 
that the granular temperature is locally determined by 
the balance between the viscous heating and the energy 
dissipation. Equation p6|) with Eqs. (fT4|) and (fT5|) gives 



T 



2 hi'') -2 
mcr , , , 7 . 



(17) 



Substituting Eq. ^ into Eqs. ^ and we get 



(18) 
(19) 



which are exactly what we have anticipated from the Bag- 
nold scaling Eqs. ((T|) and 0. 

The above derivation of the Bagnold scaling by the 
kinetic theory gives the definite expression for Eq. ([3]), 



(20) 



N{iy, T) 



(13) 



as a function of the packing fraction u. This is plotted in 
Fig.[T]by symbols connected by lines, along with the sim- 
ulation data. One can see clear discrepancy between the 
theory and the simulation especially in the higher den- 
sity region. The kinetic theory gives increasing functions 
S/N of V, which means that the flow down steeper slope 
is denser. 



III. SIMULATIONS 

In this section, we compare the expressions Eqs. (jl3p - 
(fTS)) with the simulation results of simple shear flow of 
inelastic hard spheres. 



A. Simulation setup 

The simulation is done under the constant volume con- 
dition with a uniform shear in a rectangular box of the 
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size Lx X Ly X Lz- The shear is appHed by the Lees- 
Edwards shearing periodic boundary conditions in the z 
direction [2^; The periodic boundary condition is em- 
ployed in the x and y directions. We employ the event 
driven method, using the fast algorithm developed by 
Isobe [13]. 

A steady shear flow with the mean velocity u{r) = 
(72,0,0) is prepared as follows. First, a random config- 
uration is prepared by the compressingprocedure pro- 
posed by Lubachevsky and Stillinger {2^ in the elastic 
system without shear under the periodic boundary con- 
dition. Secondly, the initial shear flow is constructed 
from the above random configuration by giving the ini- 
tial mean velocity u{r) = (72,0, 0) and setting the ini- 
tial temperature T « 100mcr^7^. Lastly, the steady 
shear flow of the inelastic system is obtained by relaxing 
the initial flow under the Lees-Edwards shearing periodic 
boundary condition [2^ . 

With the present parameter and system size, the final 
steady state is the simple shear flow with uniform packing 
fraction z/ = i/q and mean velocity u = (72,0,0) [30| . 
All the following data are taken in the steady state, and 
averaged over the space and time (typically over 10,000 
collisions per particle) unless otherwise noted. 

In the following, all the quantities are given in the di- 
mensionlcss form with the unit mass m, the unit length 
a, and the unit time 7^^. Most of the data are from 
the simulations with system size — 20, Ly = 10, 
and L2 = 40. Several simulations has been done with 
Lx — Ly — Lz — 40 to check the system size effect. We 
measure the temperature T, the normal stress N, the 
shear stress S, and the energy dissipation rate F for vari- 
ous values of the packing fraction v. These are compared 
with Eqs. p^ - p3|) from the kinetic theory. 



B. Simulation results 



1. The radial distribution function 
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FIG. 2: (color online) (a)The radial distribution functions for 
v — 0.58 with Bp = 0.7 (dashed line), 0.92 (dotted line), and 
0.98 (solid line). The peak values of contact are 136 (out of 
range), 30.6 (•), and 17.4 (o) for ep=0.70, 0.92, and 0.98, 
respectively. The theoretical value at contact, go, t (0.58), is 
shown by an arrow, (b) Plot of go,ni{i^) versus the packing 
fraction u for Cp = 0.7 (■), 0.92 (•), and 0.98 (o). go,T{i^) is 
shown by a solid line. 



As for the inelastic hard spheres under shear, a gen- 
erally accepted form of (?o(^) does not exist, but it has 
been found in several simulations that go{i') is larger for 
stronger inelasticity (l5l.[l6t. Figure[2fa) shows the radial 
distribution g{r) averaged over the all directions obtained 
from our shear flow simulation with the packing fraction 
v — 0.58 for various values of Cp. The spatial mesh to 
measure g{r) was taken as 0.001, and the peak values of 
g{r) around r = 1 (at the distance of the particle diam- 
eter) are marked by symbols for Cp = 0.98 and 0.92. We 
can see that the peak value strongly increases for smaller 
Cp, and can be much larger than the value from Eq. (j22p 
(5o,t(0.58) shown by an arrow). It is quite difhcult to 
evaluate the precise value of go{i') from this direct mea- 
surement of g(r) because of the strong increase of g{r) in 
the limit of r -|-1. 

The way we determine go{i') from the simulation is 
through the expression of the collision frequency loq (bsI . 
[36j from the kinetic theory [ll|, [s^l , 



For the constitutive relations with Table [H we need 
to know the radial distribution function at the particle 
diameter, g(){v), as a function of the packing fraction. 
For elastic hard spheres [cp = 1) in equilibrium, the well 
known expression of gQ{v) is the Carnahan-Starling for- 
mula ii 



Mv, T) - 2^g^{v)Vfvn-^'^ - ^c*(ep) j , (23) 

where c*(ep) is given in Table |T1 By measuring c^o and T 
for each v from the simulation, we can evaluate 



go,Cs{v) = 



1 - iy/2 
(1^ 



(21) 



for < v < Vf, where Vf is the freezing packing fraction 
and Vf « 0.49 Torquato [13] proposed the formula 
that include the higher packing fraction up to the random 
closed packing fraction Vc ~ 0.64 as 



6'o,t(i^) 



1 .9o,cs(i^) 
I 5o,cs(j^/)(j^/ 



for 
for 



< ly < Vf, 
(22) 



24(l-c*(ep)/32)VTi^ 



(24) 



go,ra{v]T,LOQ) is plotted versus v for various values of Cp 
in Fig.[2)Jb), where 50, t('^) in Eq. is shown by a solid 
line for reference. (7o,m('^; T, wq) shows stronger increase 
upon increasing the packing fraction v as Cp gets smaller; 
by comparing it with Fig. [Ifa), we see that this indirect 
estimate gives an reasonable Bp dependence of goii^)- In 
the following, we use go^tn{h';T,uJo) as go{v) in Table U 
unless otherwise noted. 



5 



(a) 
120 

40 



ep=0.98 s™"lation 


1 1 1 ' 
o i' : 


theory 


-X- 1 7- 


« simulation 


• ' ? / 






simulation 




_ep=0.70,Heo,y 


X/' 








-'fi'^ 1 



(b) 



600 

400 
z 

200 



1 1 1 1 1 1 
ep=0.98 simulation 


^1 


1 1 . 


theory 






« simulation 


• 




-0 =0 92 
P theory 




7 - 


„ _„ simulation 


■ 




ep=0-70 theory 





















0.50 



0.60 



0.50 



0.55 



0.60 



0.4 



0.5 

V 





1 1 1 1 1 1 1 


1 


- 8 


_ 1 1 1 1 1 1 1 1 1 ^ 

: ep=0.98 r^^Xz 




- 6 


- ^ ' — 

'~"\ ~\ 1 1 1 1 1 1 1 ~ 




_ 0.4 0.5 0.6_ 
ep=d.92* 


- ep=0.70 

* 





0.6 



FIG. 3: (color online) The energy dissipation rate F (a) and 
the normal stress A'^ (b). The simulation data for = 0.98 
(o), 0.92 (•), and 0.70 (■) are compared with the values from 
the kinetic theory (F(;/, T) and N{i',T)) shown by symbols 
connected by dashed lines for Cp = 0.98 (x), 0.92 (+), and 
0.70 (*). r(i/, Tcoii) and N{v,Tcou) with 5o,m(j^; Tcoii, t^o) are 
denoted by the solid lines, which agree with the simulation 
data (see text). 



2. The energy dissipation rate and the normal stress as 
functions of the packing fraction 

In Fig. [3l the energy dissipation rate F (a) and the 
normal stress N (b) are shown for various values of the 
packing fraction v and the restitution coefficient Bp. For 
the normal stress, we find in the simulation that Na de- 
pends on the direction a, but the differences among them 
are at most 10% in the plotted region and are not signif- 
icant compared to the difference from the kinetic theory 
that we will study in the following. Thus, here we plot 
the average N = (N^ + Ny + N^)/3. 

The value from the kinetic theory are shown in 
Fig. [Hl^a) and (b) by symbols connected by dashed lines. 
We see in Fig. (SJa) that the energy dissipation rate is 
overestimated by the theory in the dense region, and the 
disagreement is larger for smaller Cp. The normal stress 
in Fig. mb) also shows a similar tendency, although the 
relative disagreements are smaller than those in the en- 
ergy dissipation rate F. 



3. The pre-collisional velocity correlation effects and the 
collisional temperature 

a. The energy dissipation. We first focus on the dis- 
crepancy in the energy dissipation rate F. From the col- 
lision rule Eq. ([5]), the energy dissipated per collision is 
given by 



(1 _ e^]-c^ -. 



(25) 



where Cn,ij = [(c^ — Cj) ■ n.^] is the relative normal veloc- 
ity of colliding particles just before the collision. Thus, 
F is given by 



F =< l^E.j > 
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FIG. 4: (color online) The temperature T and "the collisional 
temperature" Tcoii vs. the packing fraction v. T and Tcoii 
are denoted by ■ and * with the dashed lines for Cp = 0.7, 
respectively, and by • and + with the dashed lines for Cp = 
0.92. The inset shows T and Tcoii for Cp = 0.98 represented 
by o and X with the dashed lines. 



Here, < A >coii denotes the average of a quantity A 
over all collisions; if the value of A is Ak at the k-th 
collision, < A >coii= Y.k=i ^fc/A^coU, where A'^coU is the 
total number of collisions. Note that Eq. (|26|) is the exact 
expression for F. 

On the other hand, the expression (jlSp from the kinetic 
theory with Eq. gives 



F(j.,T) = (l-e2).T 



1 + 3c*(ep)/32 



1 - c*(ep)/32 



^ncooiiy,T). (27) 



To interpret this expression, let us consider the random 
collision of particles whose velocity fluctuation is given 
by the Maxwellian. In this case, j < c^ ^j >coii= T, then 
Eq. ^ gives 



1 



r = (1 - ep)T-nuJo 



(28) 



(26) 



The difference between this and Eq. (p7|) comes from 
the deviation of the velocity distribution from the 
Maxwellian, but the difference is found to be small in the 
parameter region studied in the present paper. There- 
fore, from the comparison of the exact expression (j26p 
with the kinetic theory expression Eq. (|27p. we conclude 
that the deviation found in Fig. [3^a) comes from the fact 
that i < cl ,j >coii< T. 

b. The collisional temperature. To confirm this 
idea, we define "the colhsional temperature" Tcoii 
^n,ij >coii /4. Figure m shows Tcoii and T as functions of 
ly. One can see that Tcou is substantially smaller than T 
for > 0.5 as is concluded above. 

To demonstrate that the discrepancy is actually re- 
solved by Tcoii, we plot F(j/,r) of Eq. ^ with 
go.m{i^,T,ujo) of (Ell) in 73(1/) replacing T by Tcoii (the 
solid lines in Fig[3fa)); This is equivalent to replace T 
in Eq. ((?7)) with Tcoii and use the measured value of the 
collision frequency for ujq. The agreement is quite good. 

c. Normal stress. Now, we consider the effect of 
TcoU < r on the normal stress N. The value of Cn,ij 
should also play an important role in the collisional com- 
ponent of the normal stress A'coii, because Cn,ij directly 
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FIG. 5: (color online) The distribution of relative normal ve- 
locity just before the collision c„. Paii(cn) and 

interval J 

are compared along with Prand(c„; Tcoii) and P,and(c„; T) for 
Cp = 0.92 with u = 0.40 (a) and f = 0.58 (b). We see that the 
difference between Paii(cn) (o) and Pintcrvai (cn) (□) is small 
for v — 0.40, but for = 0.58, the Paii(c„) has sharper distri- 
bution. The solid lines show Prand (cn ; rcoii) and the dashed 
lines show Prand(c„;r). See text for details. 




FIG. 6: (color online) The pre-collisional relative normal ve- 
locity distributions Paii(c„) (o) and Pmtervai (cn) (□) com- 
pared along with Prand(c„; Tcoii) (solid line) and Prand(cn;r) 
(dashed line), respectively, for ep = 0.70 with u = 0.58. 



determines the momentum transfer from the particle i to 
the particle j through a coUision: Apj = {c* — Cj) — [{1 + 



-Api. Thus, we expect that A^coii is ap- 



proximately proportional to < |Ap| > nujQ oc ^/T^nuo. 
In addition, the coUisional part A^coii is dominant in the 
dense region. 

In Fig. Mh), iV(z^,Tcou) of Eq. ^ is plotted by the 
solid lines, where 5o,in(i^; Tcoii, i^o) is used as go(i^)- We 
see that the solid lines show reasonably good agreement 
with the data for the whole density region. 



collide with same colliding partners inelastically many 
times within a short period of time. Under the shear, 
however, this correlation will be lost when they are forced 
to pass each other and collide with new partners. The 
typical time scale that a pair of particles pass each other 
is the unit time, i.e., 7~^. This argument explains the 
smaller Tcou in the denser region, because particles collide 
more frequently with same partners before they move far 
apart |38[. 

This argument tells that the collision does not have 
memories of the previous collisions earlier than the unit 
time 7~^. To confirm this, we compare the following two 
distributions of the pre-collisional velocity: (i)i-'aii(cn), 



which is the distribution of i 



for all collisions between 



all pairs of particles, and (ii) Pintcrvai (cn), which is the 
distribution of Cn.ij of the collisions whose colliding pairs 
of particles did not collide with each other during the last 
unit time 7~^. If the velocity correlation mainly comes 
from the multiple collision with same partners within the 
unit time scale, then -Pintcrvai(cn) should have the width 
determined not by TcoU but by the average temperature 
T. 

The results are shown in Fig. [5] for = 0.92 with ly = 
0.40 (a) and 0.58 (b), where Paii(c„) is denoted by o and 
-Pintorvai(c„) is denoted by □. We see that Pintcrvai(c„) is 
wider than Pan(cn) for the denser case (b). 

If the particles with the Maxwellian velocity distribu- 
tion with temperature T collide among themselves ran- 
domly, the distribution of c„ is given by 



-Prand(c„;f ) = ^ ^Xp 



4T 



(29) 



In Fig. El Prand (c„;T) and Prand(c„; Tcoii) are shown 
by the dashed and solid lines, respectively. They are 
indistinguishable for v = 0.40 in Fig. ^a), but show 
clear difference for — 0.58 in Fig. [H^b). We find 

that Prand (Cn;r) fitS Pintcrval (Cn) , and Prand (c„ ; Tcoll) fitS 

Pan(c„), which further confirms that colliding partners 
are correlated in the way characterized by Pcoii- 

For smaller e^, the shape of the distributions devi- 
ates from Eq. based on the random collision. Fig- 
ure [S] shows Paii(c„) and Pintcrvai(cn) Compared along 
with Prand (c„; Tcoll) and Prand (c„;T), respectively, for 
Cp — 0.70 with V = 0.58. Paii(c„) has sharper distribu- 
tion than Pintervai (cn), but neither of them fit well with 
-Prand (c„; Tcoll) nor Prand(c„;T). This suggests stronger 
correlation than the case of Cp ~ 0.92. 



4- Origin of the pre-collisional velocity correlation 

One of the possible origins of the pre-collisional veloc- 
ity correlation that makes Tcoii < T is the inelasticity, 
which makes the relative normal velocity smaller upon 
collision. In this subsection, we examine how the pre- 
collisional velocity correlation develops in the shear flow. 

It is expected that the correlation grows when particles 



5. The spatial correlation in the velocity fluctuation 

To understand the the velocity correlations in more 
detail, we study the spatial velocity correlation function 
defined as 

< T,^.J [ca,iC0jA{R - |(n - rj) ■ e^l)] > 



<E„-[A(i?-|(r.-r,)-e^|)] > 



(30) 
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FIG. 7: (color online) The spatial velocity correlation func- 
tions for Bp — 0.92 and f = 0.55. (a) The correlations in the 
X direction C^^p{R). One sees that the longitudinal compo- 
nent Ci^l{R) (solid line) has the larger amplitude than others 
(shown by dashed lines), (b) The longitudinal velocity corre- 
lations Ca,a{R) btT a = X (solid line), y (dashed line), and z 
(dotted line). 



where a, /3, and 7 take x^y or z, Ca^i = (c^^i — Mq.), 
< • • • > denotes the time average, and represents the 
unit vector in the 7 direction. A{R — \r ■ e^l) is one 
when \R ~ \r ■ BjW < 0.05 and \r • e^/| < 0.1 for 7' ^ 7, 

and it is zero otherwise. We calculated C^^}{R) for the 
system with Cp — 0.92, both for the small system with 
Lx = 20, Ly = 10, Lz — 40 and the large system with 
Lx 40, Ly = 40, Lz = 40. We find that the correlation 
extends over the whole system in the case of the small 
system, but it goes to zero for the large system. In the 
following, we present the spatial correlation measured in 
the large system, but we confirmed that the hydrody- 
namic quantities presented in the previous subsections 
did not show any differences. 

In Fig. [Zl^a) , the various components of the correlation 
in the x-direction C^^^(i?) are shown. We find that the 

longitudinal correlation in the x-direction, ci'^}:{R), has 
larger amplitude than other components; this tendency 
is also found in the y- and z-direction (data not shown). 
The longitudinal correlation at the particle diameter dis- 
tance (i? = 1) is positive, which is consistent with the 
fact that Tcou < T. It is evident that the correlation 
shows an oscillation, whose wavelength is order of the 
particle diameter, which will be discussed in section HVl 
The longitudinal components in x, y and z directions 
are shown in Fig. \T[h). All of them show oscillations in 
the particle diameter scale. We also found that the lon- 
gitudinal correlation shows larger amplitude for smaller 
restitution coefhcient Bp and/or larger packing fraction v 
(data not shown). 

6. The packing fraction dependence of the shear stress 

We find that the shear stress S shows more complicated 
packing fraction v dependence than those of the energy 
dissipation rate F and the normal stress N. In FigjSl 



FIG. 8: (color online) The shear stress S vs. the packing 
fraction u. The simulation data are plotted for = 0.98 
(o), 0.92 (•), and 0.70 (■). The kinetic theory constitutive 
relations Siu, T) is shown by the symbols connected by dashed 
lines (x for = 0.98, -f for 0.92, and * for = 0.70). 



the simulation data of the shear stress S are denoted by 
symbols, and S{v^ T) from the kinetic theory (Eq. (fT4|) 
with Table [l| are denoted by symbols with the dashed 
lines. We find that, for Cp = 0.98, the shear stress is 
underestimated by the theory, while for gp = 0.92 and 
0.70, the shear stress is overestimated. 

Actually, in the case of the elastic (cp — 1) hard sphere 
system, the Enskog theory is known to underestimate 
the shear viscosity in dense region [2^, [40[ , and this ten- 
dency is seen in the result for ep = 0.98. The results for 
Cp = 0.70 shows that the inelasticity reduces the shear 
stress to the value smaller than the one expected from 
the kinetic theory, but we do not understand the reason 
of this reduction yet. Rather good agreement in between 
for the case of er, = 0.92 seems to be accidental. 



IV. DISCUSSION AND SUMMARY 

A. The shear stress and the anisotropic correlation 

In contrast to the energy dissipation rate F and the 
normal stress iV, the discrepancy in the shear stress 
S cannot be understood just by the pre-collisional ve- 
locity distribution averaged over all directions, but the 
anisotropy of the pre-collisional correlations in both the 
velocity and the position should by important in the 
shear stress. These anisotropics are not taken into 
account in the kinetic theory employed in the preset 
analysis. In fact, for the soft-sphere system in two- 
dimensional, sheared flow, it has been found that the 
contact force distribution strongly depends on direction 
[4l| . Our preliminary results also show a similar direc- 
tion dependence in the coUisional momentum transfer per 
unit time. The detailed analysis is left for future studies. 
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B. The packing fraction dependence of the ratio of 
the shear stress to the normal stress 

As we have seen in Fig. [l] S/N in the simulation is a 
decreasing function of the packing fraction v for larger 
packing fraction v and/or smaller restitution coefficient 
Cp, while Eq. ((20|) . from the kinetic theory, S/N always 
increases with v. 

Kumaran argued that the particle roughness is neces- 
sary for S/N to have a decreasing part upon increasing 
1/ in the dense region [Tsj . However, even for the smooth 
particles, the present simulations show that S/N has a 
decreasing part in the dense region for the inelastic hard 
sphere system, although the particle roughness may well 
amplify the decreasing part of S/N. 

The present authors have suggested [61] that the origin 
that leads the kinetic theory to the increasing S/N on i/ 
even for the denser region is that fs^v) in the energy dis- 
sipation r of Eq. (fT5)) increases too sharply for larger v. 
In this paper in section [III B 31 we showed that the sharp 
increase in T can be weaken by using Tcoii instead of T. In 
the present treatment, however. It is not possible to ex- 
tract the i/-dependence out of T{i', Tcoii) and to compare 
it directly with fsii') because T and Tcou are determined 
by v and 7 in the steady state simulations, therefore, the 
quantity that corresponds to fsiv) in eq. (|15|) cannot be 
defined from the simulation data. 

Finally, let us comment on the fact that S/N does 
increase with v ina certain parameter range in our simple 
shear flow simulation, in contrast to the fact that the 
increasing v upon increasing S/N = tan 9 has never been 
observed in the granular flow down a slope. This suggests 
that the steady flow in this parameter region is unstable 
in the slope flow configuration. It is interesting to study 
the relation between the stability of the flow and the v 
dependence of S/N. 

C. Oscillation in the spatial velocity correlation 

As shown in Fig. [7l the spatial velocity correlation is 
found to oscillate in the scale of the particle diameter. 
Although we have not yet understood the origin of this 
oscillation, it is plausible that the oscillation comes from 



the coupling between the density correlation and the ve- 
locity correlation. Analysis on the sheared Langevin sys- 
tem suggests that the spatial velocity correlation is re- 
lated to the radial distribution function (39j . which os- 
cillates in the particle diameter scale. It is likely that 
similar coupling also exists in the granular shear flow. 

D. Summary 

We have simulated the simple shear flow of the smooth 
inelastic hard sphere system by molecular dynamics sim- 
ulations. We have found that the energy dissipation rate 
r and the normal stress N are smaller than those ex- 
pected from the kinetic theory. We have showed that the 
relative pre-collisional normal velocity of colliding pairs 
of particles, Cn,ij, is smaller than the one expected from 
random collisions, and this reduces F and N. By exam- 
ining the distributions of Cnjj for all collisions (Paii(cn)) 
and for only the flrst collisions of the new pairs during the 
last period of time (^Intcrvailcn)), we have concluded 
that the reduction of the relative velocity is caused by the 
multiple inelastic collisions during the time period 7~^. 

To understand the velocity correlation in more detail, 
we have studied the spatial velocity correlation. It has 
been found that the longitudinal components of the cor- 
relations have larger amplitude with the oscillation in the 
scale of the particle diameter. 

The shear stress S has been found to be overestimated 
for smaller Cp, but underestimated for larger Cp by the 
kinetic theory. 
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